# Nighttime lights and land cover data are in a single chunk to avoid external pointer errors when rendering
#---
# Nighttime lights
#---
# Subset nighttime lights rasters to each city
city_sf_list <- list(flagstaff_sf, homerglen_sf, beecave_sf,
scranton_sf, belvidere_sf, clanton_sf)
city_ntl_rasters <- vector("list", length = length(city_names))
# Load files for each year
for (i in 1:length(years)){
pb_download(paste("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_t",
years[i], ".tif", sep = ""), repo = "cstaebell/final-project-cstaebell")
rast_file <- paste("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_t",
years[i], ".tif", sep = "")
year_rast <- rast(rast_file)
# Loop over cities
for (j in 1:length(city_sf_list)) {
# Crop and mask to city extent
current_rast <- year_rast |>
terra::crop(city_sf_list[[j]]) |>
mask(city_sf_list[[j]])
# Create list entry during first iteration, then append in subsequent iterations
if (i == 1) {
city_ntl_rasters[[j]] <- list(current_rast)
} else {
city_ntl_rasters[[j]][[length(city_ntl_rasters[[j]]) + 1]] <- current_rast
}
}
}
# Stack city lists
city_ntl_rasters <- lapply(city_ntl_rasters, rast)
# Create separate rasters for each city
#city_ntl_var_names <- paste0(city_names, "_ntl_rast")
#for (c in 1:length(city_ntl_var_names)) {
# assign(city_ntl_var_names[c], city_ntl_rasters[[c]])
#}
# Calculate annual means
annual_ntl_means <- vector("list", length = length(city_names))
for (i in 1:length(city_ntl_rasters)) {
ntl_vals <- data.frame(values(city_ntl_rasters[[i]])) |>
summarize(across(t2014:t2024, ~ mean(.x, na.rm = TRUE))) |>
pivot_longer(cols = t2014:t2024) |>
mutate(name = years, city = city_names[i]) |>
rename(year = name, mean_ntl = value)
annual_ntl_means[[i]] <- ntl_vals
}
ntl_means_df <- list_rbind(annual_ntl_means) |>
left_join(all_annual_acs, by = join_by(city == city, year == year)) |>
mutate(percap_mean_ntl = mean_ntl / pop_est)
#---
# Land cover data
#---
# Download land cover data
pb_download("MCD12Q1.061_LC_Type1_doy2019001000000_aid0001.tif",
repo = "cstaebell/final-project-cstaebell")
pb_download("MCD12Q1.061_LC_Type1_doy2020001000000_aid0001.tif",
repo = "cstaebell/final-project-cstaebell")
land_use_rast_2019 <- rast("MCD12Q1.061_LC_Type1_doy2019001000000_aid0001.tif")
land_use_rast_2020 <- rast("MCD12Q1.061_LC_Type1_doy2020001000000_aid0001.tif")
# Subset land cover rasters to each city
land_use_rasters <- vector("list", length = length(city_sf_list))
urban_ntl_dfs <- vector("list", length = length(city_sf_list))
urban_ntl_vals <- vector("list", length = length(city_sf_list))
group_assign <- c("Dark Sky", "Dark Sky", "Dark Sky", "Comparison", "Comparison", "Comparison")
pair_assign <- c("A", "B", "C", "A", "B", "C")
for (i in 1:length(city_names)) {
# Specify which data to use (2019 for all cities but Clanton)
if (i == 6) {
land_use_rast <- land_use_rast_2020
} else {
land_use_rast <- land_use_rast_2019
}
# Crop and mask to city extent
current_rast <- land_use_rast |>
terra::crop(city_sf_list[[i]]) |>
mask(city_sf_list[[i]])
# Save to list for future use
land_use_rasters[[i]] <- current_rast
# Mask all but urban land use (land use code = 13)
urban_area <- current_rast
urban_area[urban_area != 13] <- NA
city_urban <- city_ntl_rasters[[i]] |>
mask(urban_area)
# Calculate mean NTL for urban areas
urban_ntl_dfs[[i]] <- data.frame(values(city_urban)) |>
summarize(across(t2014:t2024, ~ mean(.x, na.rm = TRUE))) |>
pivot_longer(cols = t2014:t2024) |>
mutate(name = years, city = city_names[i]) |>
rename(year = name, mean_urban_ntl = value)
# Compile values for statistical tests
urban_ntl_vals[[i]] <- data.frame(values(city_urban)) |>
pivot_longer(cols = t2014:t2024) |>
filter(!is.na(value)) |>
mutate(city = city_names[[i]], group = group_assign[i], pair = pair_assign[i])
}
# Combine all mean NTL data into single data frame
all_ntl_means_df <- ntl_means_df |>
left_join(list_rbind(urban_ntl_dfs), join_by(year == year, city == city)) |>
mutate(mean_percap_urban_ntl = mean_urban_ntl / pop_est,
pair = case_when(
city == "flagstaff" | city == "scranton" ~ "A",
city == "homerglen" | city == "belvidere" ~ "B",
city == "beecave" | city == "clanton" ~ "C"),
group = case_when(
city == "flagstaff" | city == "homerglen" | city == "beecave" ~ "Dark Sky",
city == "scranton" | city == "belvidere" | city == "clanton" ~ "Comparison")
)
# Prepare data for sample land cover plots
# Define land cover encoding
Land_Cover_Type_1 <- c(
"Water" = 0, "Evergreen Needleleaf forest" = 1, "Evergreen Broadleaf forest" = 2,
"Deciduous Needleleaf forest" = 3, "Deciduous Broadleaf forest" = 4, "Mixed forest" = 5,
"Closed shrublands" = 6, "Open shrublands" = 7, "Woody savannas" = 8, "Savannas" = 9,
"Grasslands" = 10, "Permanent wetlands" = 11, "Croplands" = 12, "Urban & built-up" = 13,
"Cropland/Natural vegetation mosaic" = 14, "Snow & ice" = 15, "Barren/Sparsely vegetated" = 16
)
lcd <- data.frame(
ID = Land_Cover_Type_1,
landcover = names(Land_Cover_Type_1),
stringsAsFactors = FALSE)
# Create data frames for plotting
flagstaff_lulc <- as.data.frame(land_use_rasters[[1]], xy = TRUE) |>
left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID)) |>
rename(lc_code = MCD12Q1.061_LC_Type1_doy2019001000000_aid0001) |>
mutate(city = "Flagstaff")
scranton_lulc <- as.data.frame(land_use_rasters[[4]], xy = TRUE) |>
left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID)) |>
rename(lc_code = MCD12Q1.061_LC_Type1_doy2019001000000_aid0001) |>
mutate(city = "Scranton")
homerglen_lulc <- as.data.frame(land_use_rasters[[2]], xy = TRUE) |>
left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID)) |>
rename(lc_code = MCD12Q1.061_LC_Type1_doy2019001000000_aid0001) |>
mutate(city = "Homer Glen")
belvidere_lulc <- as.data.frame(land_use_rasters[[5]], xy = TRUE) |>
left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID)) |>
rename(lc_code = MCD12Q1.061_LC_Type1_doy2019001000000_aid0001) |>
mutate(city = "Belvidere")
beecave_lulc <- as.data.frame(land_use_rasters[[3]], xy = TRUE) |>
left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2019001000000_aid0001 == ID)) |>
rename(lc_code = MCD12Q1.061_LC_Type1_doy2019001000000_aid0001) |>
mutate(city = "Bee Cave")
clanton_lulc <- as.data.frame(land_use_rasters[[6]], xy = TRUE) |>
left_join(lcd, by = join_by(MCD12Q1.061_LC_Type1_doy2020001000000_aid0001 == ID)) |>
rename(lc_code = MCD12Q1.061_LC_Type1_doy2020001000000_aid0001) |>
mutate(city = "Clanton")
all_lulc <- bind_rows(flagstaff_lulc, scranton_lulc, homerglen_lulc, belvidere_lulc,
beecave_lulc, clanton_lulc)
all_lulc$city <- factor(all_lulc$city, levels = c("Flagstaff", "Scranton",
"Homer Glen", "Belvidere",
"Bee Cave", "Clanton"))